Efficient ff{N^) approach to solve the Bethe-Salpeter equation for excitonic bound states 



F. Fuchs, C. Rodl, A. Schleife, and F. Bechstedt 
Institut fiir Festkorpertheorie und -optik, Friedrich-Schiller-Universitdt and European 
Theoretical Spectroscopy Facility (ETSF), Max-Wien-Platz 1, 07743 Jena, Germany 

(Dated: May 6, 2008) 

Excitonic effects in optical spectra and electron-hole pair excitations are described by solutions of the Bethe- 
Salpeter equation (BSE) that accounts for the Coulomb interaction of excited electron-hole pairs. Although 
for the computation of excitonic optical spectra in an extended frequency range efficient methods are available, 
the determination and analysis of individual exciton states still requires the diagonalization of the electron- 
hole Hamiltonian H. We present a numerically efficient approach for the calculation of exciton states with 
quadratically scaling complexity, which significantly diminishes the computational costs compared to the com- 
monly used cubically scaling direct-diagonalization schemes. The accuracy and performance of this approach 
is demonstrated by solving the BSE numerically for the Wannier-Mott two-band model in k space and the semi- 
conductors MgO and InN. For the convergence with respect to the k-point sampling a general trend is identified, 
which can be used to extrapolate converged results for the binding energies of the lowest bound states. 

PACS numbers: 71.10.Li, 71.15.Dx, 71.15.Qe, 71.35Cc, 71.70.Gm 



I. INTRODUCTION 

The Coulomb interaction of excited electrons and holes 
plays an important role for the optical properties of condensed 
matter Ul B ISD ■ Photon-induced two-particle electronic exci- 
tations are accompanied by the rearrangement of the remain- 
ing electrons in a solid, so that the individual particles are 
renormalized to quasiparticles (QPs). Their description within 
the framework of many -body perturbation theory (MBPT) has 
made substantial progress in the last three decades fl. The 
most common approach is Hedin's GW approximation 
which describes the response of the remaining electrons by a 
dynamically screened Coulomb potential W . The self-energy 
operator E of an excited particle is thereby given by E = itiGW 
with the single-particle Green's function G. Its numerical im- 
plementation [6, 7] usually yields single-particle excitation 
energies in good agreement with angle-resolved or inverse 
photoemission experiments |^ [l^ . However, the optically 
excited quasielectron-quasihole pairs show additional interac- 
tions. They are described by the so called polarization func- 
tion P, which obeys a Bethe-Salpeter equation (BSE) ^ 01 • 
Apart from the bare electron-hole exchange, which can be 
identified with crystal local-field effects (LFEs) |2], its ker- 
nel is given by the derivative — (//j)^'5E/5G which is usu- 
ally replaced by the leading term —W. It clearly represents 
the screened Coulomb attraction between quasielectrons and 
quasiholes [3]. 

Optical spectra of real materials can be calculated from first 
principles by solving the eigenvalue problem (EVP) for an 
effective two-particle Hamiltonian H corresponding to a re- 
formulated BSE. The eigensystem of H then can be used to 
obtain a spectral representation of P. Usually, one of the fol- 
lowing numerical approaches is applied to calculate P or the 
directly related frequency-dependent dielectric function e{(o): 
(i) The explicit diagonalization of H (solving the EVP di- 
rectly ifTH [T2I1 ). (ii) the iterative Hay dock method f\3\, or 
(iii) a time-evolution algorithm that is based on an initial- 
value formulation of the Fourier-transformed dielectric func- 
tion ifTill . Meanwhile, the applications extend also to rather 



complex materials and structures such as semiconductor sur- 
faces ifsl [T6ll . nanocrystals and molecules fl?', Ts*], or even 
ice and water ifioi I20I1 . The most important changes in the 
optical spectra with respect to the independent-particle limit 
concern a general redshift of the transition energies and a re- 
distribution of oscillator strength towards lower energies iffill . 

For some systems also bound states of the electron-hole 
pairs - so called excitons fl^ - form and can be observed in 
optical spectra below the QP absorption edge. Thereby, two 
fundamentally different types of excitons - the Frenkel and 
the Wannier-Mott type - are distinguished. Excitons of the 
Frenkel type emerge in systems where the gap is confined by 
rather localized electronic states, as it is found for surfaces of 
covalently bonded semiconductors [IS] , insulators with strong 
ionic bonds llT2llT3ll22ll23ll . or crystals with hydrogen-bridge 
bonds l!2i|l. For materials with a pronounced dispersion of 
the first conduction band excitons of the Wannier-Mott type 
can form below the gap. They give rise to a 
hydrogen-like spectrum of pair eigenvalues and, in contrast to 
the Frenkel type, an electron-hole distance much larger than 
the lattice constant. Typically, their binding energies are much 
smaller than those of Frenkel-like excitons. Nevertheless, the 
lowest excitonic states can be probed by a variety of spectro- 
scopic techniques. 

However, first-principles calculations, as outlined above, 
remain a challenging task especially for Wannier-Mott-like 
excitons with binding energies below 0. 1 e V. The most impor- 
tant limitation is due to the high k-point densities required for 
sampling the Brillouin zone (BZ), in order to describe the lo- 
calization of the excitonic wave function in the Fourier space 
sufficiently. The number of k points directly relates to the 
rank of the BSE Hamiltonian, causing extreme demands in 
terms of storage and CPU time for the calculation of the spec- 
tra. The latter one rapidly becomes the limiting factor, since, 
of the previously discussed approaches, only the direct diag- 
onalization of H - involving computational costs scaling like 
ff{N^) - can be employed in the computation of bound pair 
excitations. Attempts have been made to solve the problem 
by introducing an interpolation scheme for the electron-hole 
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interaction 1*28^ or by k-point samplings restricted to only a 
small part of the BZ Ii29t,i30,1 . Still, the question how the com- 
plete electronic structure influences the bound pair excitations 
and their the oscillator strengths remains open. Also the ef- 
fects of electron-hole exchange on the excitons lIsTl [32I1 have 
been studied little so far by a combination of ab initio elec- 
tronic structure calculations and the BSE treatment of the two- 
particle excitations. Another important point concerns the 
prediction of excitonic effects, e.g. binding energies and oscil- 
lator strengths, for semiconductors such as InN, whose sample 
quality does currently not allow their measurement. Also the 
influence of polytypism of a material remains an open ques- 
tion. 

In order to reduce the computational demands necessary 
for systematic studies addressing the aforementioned ques- 
tions, we combine the use of special k-point sets and an it- 
erative matrix-diagonaUzation scheme, obtaining a computa- 
tionally very efficient approach for the calculation of excitonic 
states below the absorption edge. The methodology of this ap- 
proach is described in detail in Sec. II. The computational de- 
tails of the underlying ab initio calculations are summarized 
in Sec. Ill, before in Sec. IV.A the approach is applied to and 
tested by means of the Wannier-Mott model exciton. The con- 
sequences of "true" ab initio band structures are studied for 
the semiconductors MgO and InN in Sec. IV.B and C. Finally, 
a summary and conclusions are given in Sec. V. 

n. METHODOLOGY 

A. Bethe-Salpeter equation and pair Hamiltonian 

The inclusion of excitonic effects in optical spectra requires 
to go beyond the independent-(quasi)particle approximation 
(IQA) for the polarization function P by taking into account 
the electron-hole attraction and exchange. For this purpose 
a BSE for the irreducible polarization can be derived from 
MBPT Clll, 



P = Po+Po{2v-W)P, 



(1) 



with the IQA polarization function Pq, the statically screened 
Coulomb potential W, and the bare Coulomb potential v with- 
out its long-range Fourier component G = 0. For a detailed 
derivation and a generalization for collinear spin polarization 
we like to refer to ||33|] . 

A commonly adopted basis for representing the BSE is 
given by the Bloch states ^„k(r) of the crystal problem, char- 
acterized by the band index n and a wave vector of the first 
BZ k. Neglecting QP effects for the wave functions, they are 
usually approximated by the wave functions obtained in the 
framework of density functional theory (DFT) from the solu- 
tion of the Kohn-Sham equation. 



(XC, Vxc) interaction. The latter potential is usually given by 
the local-density (LDA) or generalized-gradient approxima- 
tion (GGA). However, for materials where these approxima- 
tions fail, such as InN, other potentials derived, e.g. from the 
LDAh-U scheme 1,34.1 or even spatially non-local potentials in 
the framework of generalized Kohn-Sham schemes Issll might 
be used. Especially the latter constitute a good starting point 
for the pertubative calculation of single-particle QP eigenval- 
ues ^^(k) M. 

With restriction to materials with completely occupied and 
unoccupied bands (semiconductors and insulators), the BSE 
^ can be rewritten as EVP for an effective two-particle 
Hamiltonian// 111 111 . 



{2n) 



dk'H^/,, (k, k')AX, (k') = E^tf}^) , (3) 



with the pair eigenvalues and eigenvectors A^,(k), Q. the 
crystal volume, and H^z the volume of the BZ. The two- 
particle Hamiltonian can be divided into a diagonal and a non- 
diagonal part, 

4^;;;,(k,k') = i/,^.(k)5,v4c'4k' +//^;Mk,k'), (4) 

with 

H^,.i^)^EQ^ik)^Efik), (5a) 

/^;;,(k,k')--w^;;,|:,+2v^,;;,[^„ (sb) 

where 5^^' is used as short notation for (27i:)^£2^^5(k — k'). 
The first, diagonal, contribution to (|4|i is constructed from the 
difference of the QP eigenvalues of conduction and valence 
bands. It describes the excitation of non-interacting quasipar- 
ticles. The second part, given by the matrix elements of the 
statically screened Coulomb potential W, accounts for the at- 
traction of electrons and holes, 

Wc":'k' = //«frflfr>i(r)(p,v(r)lV(r,r')(Pvk(r')<p;k'(r') 

£2 ,4^, |k - k' + G| |k - k' + G'^^^-'^'^^)^^'^' )' 



■ GG 



(6) 



with the symmetrized inverse dielectric function e^g, and the 
Bloch integrals 
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iio Jq.0 



(7) 



Thereby, G denotes vectors of the reciprocal lattice, £2o the 
unit cell volume, and M„k(r) the cell-periodic part of the Bloch 
waves ^„k(r) = (i2)^^''^e''^'"M„k(r). The third part contains 
matrix elements of the bare Coulomb potential and models 
LFEs by means of an electron-hole exchange term ifsl lslll . 



2m 



A + V{r)+VH{r)+Vxcir) 



(pnkir) = e„(k)(p„k(r), v'^jI,^, = jj drdr' (p*y,{r)(p,^,{r)v{r - v')^,^y^{r')^;,^,{r') 



(2) 

with the Kohn-Sham eigenvalues e„ (k) and the potentials of 
electron-ion (V), Hartree (V//), and exchange and correlation 



(8) 
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B. Generalized eigenvalue problem 



C. Optical spectra 



For any practical purpose the k-continuous formulation of 
the BSE ([J) needs to be discretized ll28ll and the sum over 
the valence and conduction bands has to be truncated. In this 
work, the latter truncation is achieved by introducing a cutoff 
energy for the transition energies of non-interacting pairs. For 
practical reasons we define this BSE cutoff with respect 
to the single-particle energies without QP corrections, corre- 
sponding to the Kohn-Sham eigenvalues, e£ (k) — ev(k) < Ect- 
For the k discretization the BZ is divided into subvolumes 
Vk with £2bz = Lk^k- The k-discrete BSE Hamiltonian and 
eigenvectors are defined as averages over the subvolumes V^, 

^«k:=^ / diiHg{ii), (9) 



rjC V k ._ 

"c'v'k' • = 



1 



dq dq'Hy;,iq,q'). 



With these definitions the discretization of (HJ yields the gen- 
eralized EVP 

" E Hc'l^V^A^,^ - (E^ Hi?,k)A^„k, (1 1) 



c'v'k' 



Fortunately, it is easy to recast (fTTT i into a conventional EVP 



E ^^'v'k'^^^v'k' - iE^-H^vk)A 



D \aA 
cvk' 



(12) 



v'k' 



using the transformed eigenvectors A^^,^ = a]iA^^,^^ with 

Special care has to be taken for the singular G = G' = 
part (SC) of W appearing in the k = k' term of H^^^,^, : 



^^IK'l'k'] 



47re2£o-i(0) 



dq[ dq'^"'^' 

Vk Jvk \q 



'\2- 



(13) 



iiVkVk iVk 'JVk " |q-q' 

It can be evaluated analytically when the six-dimensional inte- 
gration is approximatively replaced by the three-dimensional 
integral /^^ c/q/ |qp over spheres 5k of the volume represented 
by each k point. The singularity correction decreases with in- 
creasing number of k points, converging to zero in the limit 
of vanishing k — k' distance. In the case of regular k-point 
meshes the singularity correction leads to a rigid shift of the 
eigenvalues to lower energies |37]. For inhomogeneously dis- 
tributed k-point sets, however, its value is not constant over 
all k points. Therefore, we do not adopt the spherical ap- 
proximation in order to avoid spurious effects on the absolute 
values and relative positions of the eigenvalues. Instead, we 
carry out the six-dimensional integral ( fTSl l numerically, tak- 
ing into account the actual size and shape of Vk- Since the 
numerical convergence of these integrals is utterly slow but 
perfectly linear when plotted over the sampling-point distance 
along one direction, we extrapolate the value of the singularity 
correction from the results of two calculations using 10*^ and 
20^ sampling points inside Vk. Still, the integration of (fT3T l is 
rather time-consuming and therefore only sensible in case of a 
limited number of different sizes and shapes of Vk, e.g. in the 
case of regular or hybrid k meshes as introduced in Sec. Ill El 



The diagonal components of the frequency-dependent 
macroscopic dielectric tensor emac{o)) can be obtained using 
the eigenvectors A^,,^ and energies from the solution of 
(O, 



An e^h^ 



f'j 



1 



£2 mo ^E^ ^^^^E^-^him + iyY 



(14) 



where the oscillator strengths /j^ are given by 
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^ {(pck\Pj\(pvk) ^ TA* 



(15) 



and {(pck \pj \ (Pyk) denote the matrix elements of the momen- 
tum operator in the Cartesian direction j. 

For the calculation of the dielectric function (DF), Schmidt 
et al. ifTill devised an efficient approach, scaling like i&{N^) 
with the dimension of the excitonic Hamiltonian. This ap- 
proach avoids the direct diagonalization of ft by reformulat- 
ing (fl4l l as an initial-value problem and solving that. It, how- 
ever, comes at the price of losing the explicit knowledge of 
the exciton excitation energies E^ and wave-function coef- 
ficients A^,|j. Further, a certain broadening 7 of the optical 
transitions is mandatory for the numerical convergence of the 
algorithm. For large and complex systems, where can reach 
up to 10^ or more, the computational costs for a full diagonal- 
ization of H become prohibitive due to their 0'{N^) scaling. 
This renders the time evolution the only feasible approach for 
the calculation of the DF. Another method, also aiming at the 
DF only, was proposed by Benedict et al. [13]. It is, however, 
hard to compare in terms of computational efficiency due to 
its implementation in real space. 



D. An 



I algorithm for excitons 



Often, not the DF itself but isolated exciton levels, espe- 
cially bound states with pair excitation energies E^ below the 
lowest QP gap, are of interest. Due to the limited spectral 
resolution and the possibility of vanishing oscillator strength, 
such excitations cannot be studied with the time-evolution 
method. The same holds for the analysis of excitonic wave 
functions, which requires the knowledge of the respective 
eigenvalues and eigenvectors. In those cases Eq. (O needs 
to be solved directly, raising again the problem of excessive 
computational costs for complex systems. For Wannier-Mott- 
like excitonic states near the absorption edge, usually a fairly 
small BSE cutoff can be used in order to reduce the dimen- 
sion of H |38]. On the other hand, a large number of k points 
is required in the computations to ensure a sufficiently large 
crystal volume to accommodate the exciton, which can extend 
over several thousand unit cells according to exciton localiza- 
tion radii in the range of 10 — 100 A. 

For excitons, numerical convergence is typically obtained 
for dimensions of H somewhere in the range 
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FIG. 1: (Color online) CPU time on a single Opteron 275 core 
needed for the full diagonalization (CHEEV, black filled circles) or 
the computation of the 15 lowest eigenvalues of randomly filled Her- 
mitian matrices, using the LAPACK library routine CHEEVX (black 
unfilled circles) or the CG+SR algorithm (green squares). The solid 
lines represent cubic (black) and quadratic fits (green) to the different 
data sets. 



Direct matrix diagonalization schemes can be used at the 
lower end of the aforementioned range with computational 
times of a few CPU hours. However, due to the rapid i?{N^) 
increase of computational costs, already problems of mid- 
range size are computationally too expensive for systematic 
studies. 

Given the combination of a large-rank Hamiltonian and the 
interest in only a few of the smallest eigenvalues and corre- 
sponding eigenvectors renders the problem strongly reminis- 
cent of the Kohn-Sham problem in DFT. Of course, this refers 
only to the dimensions of the problems, because in fact the 
BSE-EVP is less complicated, since no self-consistent update 
of the Hamiltonian is required. Using this observation, we 
propose to employ iterative minimization techniques similar 
to those used in DFT for the study of distinct exciton levels. 
These are based on the minimization of the Ritz functional. 



(16) 



with the trial vector | ^^). Obviously, an unconstrained mini- 
mization yields 



(17) 



the lowest (first) eigenvalue and the corresponding eigenvec- 
tor of H. For any higher eigenvalue A > 1 the minimiza- 
tion needs to be constrained to the orthogonal complement of 
span(A' , . . . ,A^^' ). For the actual calculation we follow the 
scheme devised by Kalkreuter and Simma |39]. There, con- 
secutive conjugate gradient (CG) steps are performed together 
with intermediate diagonalizations of H in the low dimen- 
sional subspace span(0^, . . . , 0"'), so called subspace rota- 
tions (SR). The dimension ris of the subspace, thereby, roughly 
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FIG. 2: (Color online) The left panel illustrates the construction 
scheme for the hybrid k-point meshes for a two-dimensional BZ. 
The k points of the coarse mesh are indicated by unfilled circles, 
those of the refined mesh by red dots. Further, the boundaries of 
the BZ and its inner part are indicated. Following the notation in- 
troduced in the text, the mesh shown here would be referred to as: 
hybrid 8:3: 21.33. The right panel shows the localization of the 
wave function for an Is Wannier-Mott exciton (with the parameters 
as in Sec. lIV Al l in reciprocal space. 



corresponds to the number of desired eigenvalues (typically 
Hj « 20). We have chosen a CG-based algorithm for the sake 
of simpUcity and due to its robustness. Other algorithms, how- 
ever, can be used as well and might also improve performance. 

As akeady stressed above, in practice only a limited num- 
ber of excitonic levels are of interest for the study of exci- 
tons. Further, this number does usually not depend on A^, the 
dimension of H. Therefore, under the assumption that the 
number of total CG steps is independent of A^, the theoretical 
operation-count scaling of the CGh-SR algorithm is found to 
be ^(a1. The reason is that it involves only matrix-vector 
and vector- vector operations of dimension N and diagonaliza- 
tions of fixed-size nj x matrices. In general, the assumption 
of a constant number of total CG steps, is of course, not true, 
since it is determined by the convergence of the algorithm. 
However, it is found that the convergence depends little on 
the matrix dimension, but more on its conditioning. Figure [T] 
confirms that the CGh-SR method scales basically like ff{N^). 
The crossover point is found below N =\ 000 when competing 
with full diagonalization, or at 3000 when comparing against 
the LAPACK |40] CHEEVX routine for calculating selected 
eigenvalues and vectors of randomly filled Hermitian matri- 
ces. 

Moreover, the CGh-SR algorithm can be easily parallelized. 
MPI parallelization can help to meet the high memory de- 
mands for solving the BSE by utilizing distributed memory. 
However, also very reasonable speedups are obtained by us- 
ing multiple processors. 



E. Hybrid k-point meshes 

Even though the algorithm suggested in the preceding sec- 
tion keeps the computational costs for the calculation of exci- 
tons tractable, one drawback of the EVP remains - its extreme 
memory demands. Using very dense k-point meshes the BSE 
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Hamiltonian matrices (|4|l easily exceed 50 GB storage. Since 
the storage requirements scale quadratic ally with the number 
of k points, also the use of refined regular k-point meshes 
rapidly approaches the limits of today's supercomputers. 

In order to reduce the memory demands, but also the com- 
puter time for setting up H, one may benefit from noting that 
the wave functions of Wannier-Mott like excitons are well lo- 
calized in k space (cf. Fig. |2]i. Therefore, it is reasonable to 
refine the k meshes only in the center of the BZ 1 12>.29.1. Such 
a refinement, however, comes at the price of having to intro- 
duce varying weights (see Eqs. dgjl-lfTTTl). turning (O into the 
generalized EVP ( fT2] i. 

The hybrid k-point meshes used in this work are con- 
structed in the following manner (see also Fig. |2|i. Firstly, 
the whole BZ is sampled using a comparably coarse mesh 
of Monkhorst-Pack k points [41J. Secondly, an inner region, 
fenced in by k points of the coarse mesh, is defined and the 
points inside this region including those on the borders are 
removed. In a third step, this inner region is filled with a 
fine-grained regular k mesh, which must be constructed ac- 
cordingly to include k points on the borders of the inner re- 
gion. Finally, the latter k points are centered inside their as- 
sociated k-space volume. In some cases it is meaningful to 
further refine the k-point mesh in the inner region. This can 
be achieved by repeating steps two and three with respect to 
the inner mesh. Hereafter, meshes constructed accordingly by 
double repetition of the latter steps are referred to as double- 
hybrid ones. 

The following notation is used to characterize the afore- 
described hybrid meshes (see also Fig.|2]and its caption). The 
first set of numbers indicates the k points in the outer (coarsely 
sampled) part of the BZ along the directions of the reciprocal 
basis. The second set indicates the size of the central part of 
the BZ (with refined sampling) with respect to the k-point dis- 
tance of the outer mesh, and the third set of numbers provides 
a measure for the sampling density inside the central part of 
the BZ. For better comparison the latter set corresponds to 
the number of regular k points in each direction necessary for 
sampling the whole BZ at the density achieved by the refined 
sampling. Double-hybrid meshes are characterized addition- 
ally by the size and sampling information of the second re- 
finement region, encoded in the same manner as for the single 
hybrid meshes described before. To shorten the notation in 
case of isotropic sampling, only one number is given per set 
(see also the caption of Fig.|2]i. 

One should note that the use of a hybrid mesh requires the 
knowledge of the dielectric screening 6^^^, (q) at all vectors 
q = k — k' for the calculation of W^"'/[^/ in This also in- 
volves q points that are not part of the hybrid k-point mesh, 
which is no concern in case of a constant or model screen- 
ing as used throughout this work. However, if the screening 
is calculated in random-phase approximation (RPA), this con- 
stitutes an additional complication, that might be overcome 
by interpolating the screening for q points not included in the 
hybrid mesh. 



III. COMPUTATIONAL DETAILS 

The results presented for MgO and InN in Sec. |IV] B and 
C are based on the results of ab initio DFT calculations, car- 
ried out using the Vienna A/? initio Simulation Package (VASP, 
142, 14311 ). The projector-augmented wave (PAW) method 
1 4^ [4511 is used to model the electron-ion interaction. Plane 
waves up to a kinetic energy of 400 eV are included in the ba- 
sis set of the electronic wave functions. Quasiparticle band 
structures are calculated for InN and MgO as a gauge by 
means of the HSEOSh-GqWo approach with computational de- 
tails as described in detail in Refs. [36l|43l Within the explicit 
calculations the single-particle eigenvalues and wave func- 
tions are computed using exchange and correlation potentials 
according to the semi-empirical LDAh-U (InN) or the GGA 
(MgO) scheme together with a scissors operator to correct for 
the gap underestimation. The BSE Hamiltonian is set up ac- 
cording to (|4]i, whereby, the inverse dielectric function enter- 
ing the screened interaction term W"/^/ is approximated by 
a G-diagonal model function ll46ll . The latter parametrically 
depends on the static electronic dielectric constant, which is 
taken from RPA calculations or experiment. 



IV. RESULTS 

In the following we present results obtained within the 
afore-introduced CGh-SR approach for three applications. 
Firstly, we numerically study the two-band Wannier-Mott ex- 
citon in k space to demonstrate the applicability of the CGh-SR 
approach and to extract some general trends. Further, the for- 
mation of excitons is studied for the large-gap material MgO 
and the narrow-gap semiconductor InN. Due to their band 
structures, the excitons occurring in both materials are ex- 
pected to show Wannier-Mott-Uke character. 



A. The Wannier-Mott model 

The probably most simplified model of a semiconduc- 
tor band structure is the two-band model of two opposed 
parabolic bands, separated by the fundamental gap Eg. Within 
this model many fundamental properties can be calculated an- 
alytically, including the formation of a series of excitons be- 
low the absorption edge as found by Wannier 147] and Mott 
fisll . The analytic solution in three dimensions VM uses the 
fact that the BSE Hamiltonian (01 for the two-band model. 



/i^kp 



4k' 



neo|k-k'|2 



(18) 

resembles that of a hydrogen-like atom with the known solu- 
tion of a Rydberg series for the discrete part of the spectrum. 



Enlm Eg R(;x/n , 



(19) 



where the quantum numbers n, I, and m of the hydrogen 
atom now label the excitonic states. Henceforth, R^x = 
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FIG. 3: (Color online) Convergence of the binding energy of the Is 
exciton with the cutoff energy Shown are the results for two reg- 
ular (squares) and five hybrid k-point meshes (circles and crosses). 
The hybrid meshes feature a common k density in the inner part, 
equivalent to that of the regular 80"' k-point mesh. The y axis refers 
to the difference of the numerically calculated binding energy and 
the analytical result of 283.45 meV. Solid lines indicate the results of 
linear fits inside the region 0.08-0.2 eV^'. Dashed lines correspond 
to extrapolations based on these fits. 



FIG. 4: (Color online) Convergence of the Is exciton binding en- 
ergy with the k-point sampling in the inner zone for three different 
sizes of the densely sampled inner zone and a fixed sampling density 
in the outer region of the BZ. The y axis refers to the difference of 
the numerically calculated binding energy and the analytical result of 
283.45 meV. The solid lines without symbols correspond to the sin- 
gularity correction in the three-dimensional spherical approximation 
( IjtIi . light blue) and the results from the six-dimensional integration 
(cf. l ll3t . dark blue) as used in the present calculations. 



7?oo/i/(moeo) denotes the excitonic Rydberg constant and Eq 
determines the effective static and spatially constant screen- 
ing of the Coulomb potential. The reduced mass jJ. — (m*^' -f 
m*r^)^^ follows from the effective masses of the valence and 
conduction band. 

For the numerical solution the effective masses are cho- 
sen to be 1.0 mo and 0.5 otq, such that /i = 1/3 mq. For the 
band gap a value of 3.0 eV is assumed. The calculations are 
performed in a simple cubic k-space volume extending over 
2n/3 ' , which naturally limits the maximum BSE cutoff to 
roughly 15.5 eV. Therefore, in the following only BSE-cutoff 
energies of 15 eV or below are considered. Using a screen- 
ing constant of Eq — 4, the excitonic Rydberg - corresponding 
to the highest possible exciton binding energy - amounts to 
7?e.v = 283.45 me V. 

In order to treat ( fTSl ) numerically in k space, it is neces- 
sary to introduce a k-point sampling. Using Monkhorst-Pack 
meshes of regular k points, we first study the cutoff depen- 
dence of the binding energy of the Is {n = 1, I — Q) exciton. 
Figure [3] shows the results for two different regular meshes 
consisting of 40 x 40 x 40 and 80 x 80 x 80 points. According 
to Fig. [3] the Is binding energy increases with an increasing 
number of states included in H and a direct proportionality 
between the binding energy and {Ecui ~Eg)^^ exists. While 
the latter relation cannot be expected to transfer to the situa- 
tion of real semiconductors with a multitude of contributing 
bands, it is useful for estimating the effect of a limited BSE 
cutoff for the two-band model. Based upon linear extrapola- 
tion of the data between Ecm = 8 — 15 eV (cf. Fig. O, the Is 
binding energy at 15 eV cutoff is expected to underestimate 



the Ecut °° value by less than 5 meV, i.e. by less than 2%. 
The energy difference between the two k-point sets is much 
larger, demanding for a better k-point sampling. 

Unfortunately, the 80 x 80 x 80 k-point Hamiltonian matrix 
at a cutoff of 10 eV already requires 95 GB RAM in single 
precision storage and would require 480 GB at Ecut — 15 eV, 
due to the relation °^ N'^'^, specific for the two-band 
model. Therefore, the hybrid k-point meshes introduced in 
Sec. Ill El are adopted in the following. The k-point density 
in the central part of the BZ corresponds to 80 x 80 x 80 k 
points in the full zone for all of them, thus allowing for a direct 
comparison to the results of the regular mesh with 80^ points. 
According to Fig. [3] it is possible to increase the distance of 
adjacent k points by a factor of two in the outer zone with- 
out introducing any significant error, if the inner zone extends 
to half of the BZ along each direction. Further, it should be 
noted that this 40 : 21 : 80 mesh contains only about a fourth 
of the k points included in the regular 80^ mesh and thereby 
reduces the memory demands by a factor of 16. If the size of 
the densely sampled inner zone is reduced further, the bind- 
ing energy increases. This is somehow unfortunate, because 
it counteracts the BSE-cutoff convergence. An analog trend is 
found, when the sampling in the outer region is reduced while 
the size and sampling density of the inner zone are kept fixed. 
Moreover, Fig. [3] shows that a trade-off between the minimum 
size of the inner zone and the sampling density of the outer 
zone exists, if one allows for a fixed error with respect to the 
corresponding regular mesh. 

Figure |4] shows the results for the Is binding energy using 
k meshes with a fixed sampling density in the outer zone and 
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FIG. 5: (Color online) Convergence of the first, second, and third 
shell (n = 1 — 3) excitons with the k-point sampling for two different 
sizes (40 : 11 - black lines with filled symbols and 40 : 7 - red lines 
with unfilled symbols) of the densely sampled inner zone and a fixed 
sampling density in the outer region of the BZ. Different angular 
momentum quantum numbers / are encoded by different symbols. 
Circles correspond to excitons with s, squares to p, triangles to d 
character. The binding energy is given with respect to the analytical 
result of 283.45 /n^ meV with n = 1 - 3. 



a varying inner density and size of the inner zone. Obviously, 
the binding energy is directly proportional to the distance of 
adjacent k points, allowing for a linear extrapolation similar 
to the situation found for the convergence with the BSE cut- 
off. We would like to point out that the observed trend is not 
dominated by the singularity correction (cf. Eq. (T3[). even 
though the latter obeys the same proportionality. This can be 
seen from Fig. |4] which also includes plots of the singularity 
correction for the three-dimensional spherical approximation 
IstIi and the six-dimensional integration (fTSl ) as utilized in the 
present calculations. Taking into account that the binding en- 
ergy and not the eigenvalue itself is plotted, it can be seen 
from Fig.|4]that the singularity coiTection akeady partially an- 
ticipates the k-point convergence. 

Figure|5]shows the k-point convergence of the first, second 
and third shell excitons for two sizes of the inner zone. In the 
realm of the two-band model all eigenvalues of a fixed princi- 
pal quantum number n are expected to be degenerate. Indeed, 
this degeneracy is found in the better-converged numerical re- 
sults with a remaining error of about 0.1 meV or less. It is 
obvious from Fig.|5] that the rate of convergence might differ 
from one shell to another, but also within a shell of equal n, 
at least in the lower sampling regime. There also the linear 
behavior can be spoiled. This is readily understood, taking 
into account that the respective wave functions differ in their 
real-space localization. Therefore, since the exciton localiza- 
tion volume increases with the principal quantum number n, 
convergence with respect to the k-point sampling is achieved 
harder for higher excitonic states. The same holds for excitons 
of lower quantum number / when comparing them to others of 
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FIG. 6: (Color online) Oscillator strengths of the n = 1 — 3 excitons. 
Vertical and horizontal dashed lines indicate the analytical results for 
the exciton binding energies and oscillator strengths. Only excitons 
of angular momentum quantum number / = are found to be bright. 
All other excitons with / > are dipole-forbidden as indicated by 
vanishing oscillator strengths. 



the same shell n. Since these states become more localized in 
the reciprocal space, it is in principle possible to use hybrid 
meshes with an even smaller inner zone for computing the lat- 
ter This also becomes apparent from Fig. |5] comparing the 
different k meshes. 

The oscillator strengths of the n = 1 — 3 excitons are 
shown in Fig. |6l From the analytic description only excitons 
with s character are expected to have non-vanishing oscilla- 
tor strengths, which decay with increasing principal quantum 
number like n^^ . Indeed the numerical results follow these ex- 
pectations with only small deviations. Further, the / > exci- 
tons are found to have vanishing oscillator strengths. This is in 
good agreement with the analytical result, predicting a direct 
proportionality between /j^ and the excitonic wave function at 
vanishing electron-hole distance in the case of k-independent 
dipole-matrix elements (cf. Eq. flSi). Since hydrogenic wave 
functions are non-vanishing at r = only for I — 0, the I > 
excitons are dipole-forbidden. 



B. MgO 

Magnesium oxide (MgO) is a wide gap semiconductor with 
a gap of approximately 7.8 eV |49]. Up to very high pressures 
its equilibrium crystal structure is given by the rocksalt (rs) 
phase 1 50]. In DFT-GGA the gap is strongly underestimated 
with a value of only 4.5 eV. Using first-order perturbation the- 
ory (Go Wo) based on the DFT eigenvalues, the gap cannot be 
fully corrected, resulting in a GGAh-GoWo QP value of 6.8 eV. 
This deficiency can be mostly overcome starting from the non- 
local HSE03 XC functional. Together with perturbative QP 
coiTections this approach results in a reliable electronic struc- 
ture with a slightly underestimated gap of 7.5 eV. Other, even 
more sophisticated QP calculations, taking into account self- 
consistency lIsTl I52I1 or even vertex corrections fsS] in the 
self-energy, are found to overestimate the gap slightly by ap- 
proximately the same amount. All of these schemes, however, 
share the drawback of being computationally by far too ex- 
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FIG. 7: (Color online) MgO band structure in the GGA+A (solid 
black lines) and HSE03+GoWo (red circles) approximation. The 
scissors shift A corrects for the gap underestimation. 



pensive to act as foundation for a study of bound excitons. 
Therefore, we use the GGA electronic structure as basis for 
the following studies, simply correcting the gap by a scissors 
shift of 2.98 eV. Figure|7]demonstrates that this approximation 
compares reasonably to the results of the more sophisticated 
HSEO3+G0W0 approach. Especially the band dispersion is 
found to compare well between both approaches. 

Due to the rocksalt symmetry of the crystal lattice, MgO 
has an i-like conduction-band minimum (CBM) of Fi t-type 
and a threefold-degenerate p-like valence-band maximum 
(VBM) of y-type. Two of the three uppermost valence 
bands show almost the same dispersion in the region around 
the F point, while the third one is more dispersive. In this situ- 
ation the formation of three degenerate, i-symmetric, excitons 
is expected from k • p theory ll53ll . In the following we will la- 
bel these excitons by A, B, and C. Indeed, for well converged 
k-point samplings A, B, and C are found to be almost degen- 
erate in our numerical calculations. The remaining splitting 
between (A, B) and C amounts to 0.30 meV, corresponding to 
less than 0.1% of the total binding energy, and might be due 
to numerical errors. In full agreement with the initial expec- 
tations all three excitons are visible in any polarization direc- 
tion. For the screening in our calculations we use an electronic 
dielectric constant of £00 = 3.0 which is in between the exper- 
imental value of 2.94 |54] and the value of 3.16 obtained in 
RPA ISOi]. The influence of £00 will be discussed below. 

In the inset of Fig. [8] the convergence with respect to the 
BSE-cutoff energy is shown for the lowest three excitons. 
Clearly, the variation is not strictly linear as in the case of 
the two band model. However, at BSE-cutoff energies higher 
than 13 eV, corresponding to (£■„„ — Eg)^^ =0.117 eV^^ the 
variation flattens out. Therefore, upon the data of Fig. |8] one 
can estimate that the exciton binding energies computed at a 
BSE-cutoff energy of 13.0 eV, as used in the following, are 
accurate within about 10 meV, i.e., less than 2.4% . 

The convergence with the number of k points for MgO is 



FIG. 8: (Color online) Convergence of the MgO A, B (unfilled sym- 
bols), and C excitons (filled symbols) with respect to the k-point 
sampling and the BSE-cutoff energy (inset). The cutoff convergence 
was studied using a regular mesh ofl6xl6xl6k points, correspond- 
ing to a minimum k-point distance of 0.16 A^'. Unfilled symbols 
refer to the binding energies of the A and B exciton, while filled 
symbols indicate those of the C exciton. The values of the hybrid 
mesh 10 : 7 are extrapolated from the region of linear variation to- 
wards zero k-point spacing. 



studied in Fig. |8] As found for the two-band model in the pre- 
ceding section, we observe a linear variation with respect to 
the distance of neighboring k points. From this data an aver- 
age binding energy for the A, B, and C excitons of 429 meV 
is derived by linear extrapolation. In comparison to experi- 
mental results of about 80 meV ||49|] and 145 meV jssl the 
calculated binding energies turn out to be drastically overes- 
timated. This may be attributed to the neglect of dynamical 
screening [5^ |57|], which can influence the exciton binding 
energies significantly. 

Since the experimental binding energies are of the order of 
the longitudinal optical phonon energies Tkolo = 89 meV ifsiil . 
the lattice polarization can partially contribute to the screeiiing 
of the electron-hole attraction. It is suggested lfT8ll57il58ll59ll 
that the static electronic dielectric constant £00 has to be re- 
placed by an effective dynamical one £0, which is enlarged 
by parts of the lattice polarizability « (fiy — £„). For MgO 
the effect of lattice polarization is very pronounced, as in- 
dicated by the large difference between the static electronic 
dielectric constant £„ = 3.0 and the static dielectric constant 
£j = 9.8 [54]. This may lead to a strong reduction of the bind- 
ing between electron and hole. An estimate of the true effect 
would reguire an explicit treatment of the dynamical screen- 
ing lf56[l57ll . which is a computationally extremely demanding 
task. Therefore, we performed test calculations using an ef- 
fective screening constant of £0 — 6.0, as suggested from fits 
to the experimental exciton spectra 15811 . finding reduced ex- 
citon binding energies of in average 99.8 meV for the A, B, 
and C exciton. 

It may be confusing to note that other theoretical studies 
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FIG. 9: (Color online) Convergence of the second shell MgO exci- 
tons with respect to the k-point sampling. Filled symbols refer to 
bright excitations, while the excitons marked by unfilled symbols 
have vanishing oscillator strength. The values of the double-hybrid 
mesh 10 : 5 : 34 : 7 are extrapolated towards zero k-point spacing. 
Further, the degeneracies of the different eigenvalue complexes are 
annotated. 



ifisi l60ll also predict much lower exciton binding energies 
closer to the experimental values, even though they use the 
static electronic dielectric constant or a RPA screened poten- 
tial as well. We attribute most of the difference to the com- 
parably coarse k-point samplings of 256 ll60ll and 1000 ifisll 
k points in the full BZ used in these studies to calculate the 
dielectric function in an extended energy range. As demon- 
strated in Fig. [8] the exciton binding energies of MgO show 
a rather strong variation with the k-point sampling density, 
leading to binding energies of 282 meV (A, B) and 191 meV 
(C) at a sampling density of 1000 k points in the BZ. 

An interesting effect of the real behavior of the electron- 
hole attraction can be observed for the higher pair excitation 
states. There, the non-spherical deviations of the potential W 
from the model |k — k'|^^ dependence, are found to cause a 
splitting of the exciton levels 4-15 into five distinct levels. 
Three of the latter are threefold-degenerate, and one of them 
is found to contain the only three excitons with non-vanishing 
oscillator strength (cf. Fig. |9]l. Figure |9] shows the conver- 
gence of these states with respect to the distance of adjacent k 
points, whereby the linear variation at small k-point distances 
indicates a reasonable convergence. 

An open question is the absolute and relative influence of 
the LFEs or electron-hole exchange effects proportional to v 
in the pair Hamiltonian H ^ on the binding of the excitons. 
In order to study this effect, we setup the BSE Hamiltonian 
according to (|4]i, once including both contributions (W and v) 
and with IV or v separately. Two k-point samplings, given 
by the hybrid meshes 10 : 7 : 24 and 10:7: 20, are used to 
account also for potentially different rates of convergence be- 
tween the two contributions W and v. As expected, excitonic 
bound states cannot form in absence of the attractive screened 



Hamiltonian 


average A, B, C binding 


energy (meV) 


extrapolated 10 : 7 : 24 


10 : 7 : 20 


H = -W + 2v 


429.2 370.5 


358.1 


H=-W 


479.6 419.9 


407.3 



TABLE I: Effect of the electron-hole exchange on the average exci- 
ton binding energy of the A, B, and C excitons in MgO. 

Coulomb interaction W between electron and hole. If, how- 
ever, the electron-hole exchange v is suppressed, the bind- 
ing energies of the lowest three excitons increase by about 
58 meV. In relation to the exciton binding energy the reduc- 
tion due to the LFEs amounts to about 14% or even more if 
dynamical screening is taken into account. The splitting of A, 
B, and C is slightly affected. From TableUit can be inferred, 
that the convergence rates of W and v indeed differ. Actu- 
ally, the contribution of v is found to be well converged at the 
probed k-point sampling densities. 

C. InN 

Indium nitride (InN) is a narrow-gap semiconductor with a 
gap of approximately 0.6-0.7 eV |.6l. .62>,63i1 . crystallizing in 
the wurtzite (wz) structure under ambient conditions. Compli- 
cations for the theoretical treatment of InN arise from the fact 
that DFT calculations including the In 4d electrons as valence 
states fail to predict a finite gap, no matter whether a local 
(LDA) or semilocal (GGA) approximation is used for the XC 
functional. Instead, the calculated band structures correspond 
to that of a zero-gap semiconductor with an, in comparison to 
the experiment, inverted ordering of the Fi ^ and , Fi ,, lev- 
els |64]. It was shown recently |36, 65] that these deficiencies 
can be overcome, using either the method of optimized ef- 
fective potentials [65.1 or a hybrid XC functional like HSE03 
ll66ll . Both methods, however, are computationally too expen- 
sive for a study of excitons. 

Therefore, we retreat to the much simpler LDAh-U scheme 
ll34l l67[ l68ll . trying to obtain a reasonable approximation of 
the HSE03h-GoWo band structure in the gap region around the 
F point by tuning the intra-atomic Coulomb repulsion of the 
In 4d electrons by the U parameter. The band structure re- 
sulting from such an LDAh-U calculation is shown in Fig.fTOl 
in comparison to the QP band structure derived within the 
HSE03h-GoWo scheme. For clarity only the results fitting 
best, obtained for U=3 eV, are shown. Moreover, a scissors 
shift of 0.48 eV is used to open the gap additionally towards 
the value of 0.71 eV, as predicted by the HSE03h-GoWo cal- 
culations. Using the LDAh-U approach, the effective masses 
at the F point amount to approximately 2.2 — 2.5 mo for the 
uppermost valence bands (va,vb) of F5.,, type. The masses 
of the crystal-field split off (vc, Fi ,,) and the first conduction 
band are much smaller with values of approximately 0.03 niQ. 
The respective values obtained from the HSE03 calculation 
differ notably, amounting to 0.9 niQ for the mass of the and 
Vfi bands and 0.046 niQ for the vc band and electron mass. 
We will try to assess the influence of the modified dispersion 
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FIG. 10: (Color online) InN band structure in the LDA+U+A (solid 
black lines) and HSE03+GoWo (red circles) approximation. A 
U parameter of 3 eV is used for a good approximation of the 
HSE03+G()W() band structure in the gap region. The scissors shift 
A corrects for the remaining gap underestimation. 



at the end of this section. The crystal-field splitting amounts 
to Ac/ =11.9 meV in the present approximation of LDA+U 
(U=3 eV), underestimating slightly the experimental values 
of 19 - 24 meV |69] but also the HSE03+GoWo value of 
22.9 meV. We do not take into account the spin-orbit split- 
ting of the valence bands, which amounts to 21.3 meV ac- 
cording to HSE03 calculations ItoIi . For the static electronic 
dielectric constant entering W a value of 7.9 is used. It 
corresponds to the average of the experimental values of the 
ordinary (7.83) iItiIi and extraordinary (8.03) ll72|] polarization 
component determined from the analysis of the respective DF 
below the wz-InN gap of 0.68 eV. 

Due to the wurtzite symmetry and the neglect of spin-orbit 
interaction, the VBM of InN is given by a twofold degenerate 
state. Therefore, the two strongest bound excitons (A, B), cor- 
responding to the Wannier-Mott Is excitons, are expected to 
be twofold degenerate as well. A third Ls-like exciton (C) is 
expected to form between the first conduction and the crystal- 
field split-off vc band, at an eigenvalue differing from that of 
A and B by roughly A^/. Since optical transitions between 
va,vb and the CBM at F are dipole-allowed only in ordinary 
polarization (E ± c) and v'c-CBM transitions only in extraor- 
dinary polarization (E || c) [73,], the corresponding excitons 
are expected to have non- vanishing oscillator strengths in the 
respective directions fl^ . 

For a moment we limit ourselves to the discussion of the 
A and B excitons, which allow for a simple definition of their 
binding energies f29'l. In agreement with the initial expecta- 
tions A and B are found to be degenerate. The inset of Fig.fTTI 
shows their convergence with respect to the BSE-cutoff en- 
ergy. The overall variation of the binding energy with respect 
to the BSE cutoff is extremely weak, due to the very small 
binding energies and the strong localization of the lowest con- 
duction band around the F point (cf. Fig [TOb . According to 
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FIG. 1 1 : (Color online) Convergence of the InN A and B excitons 
with respect to the k-point sampling and the BSE-cutoff energy (in- 
set). The cutoff convergence was studied using the double hybrid 
k-point mesh at a k-point distance of 0.0036 A^' . The values of the 
double-hybrid meshes 10 x 10x6 : 1 : 90x90x54 : 1 are extrapolated 
towards zero k-point spacing. 



Fig. [TT] a BSE cutoff of 2 eV - as used in the following - 
should introduce an error of less than 0. 1 meV. 

The convergence with the number of k points is also 
demonstrated in Fig. [TT] Due to the combination of a nar- 
row gap, large screening, and a very small effective electron 
mass found for InN, extremely dense k-point meshes are re- 
quired to converge the exciton binding energies. Therefore, 
the double hybrid k-point meshes, as introduced in Sec. Ill El 
are used to study the formation of excitonic states in InN. At 
the high-sampling limit of the studied k-point meshes the lin- 
ear variation with respect to the distance of adjacent k points, 
as observed for the two-band model and MgO, is restored. 
Extrapolating the binding energies of the lowest excitons to- 
wards zero k-point distance, reveals binding energies of about 
5.0 meV for the A and B exciton. The absolute difference 
between the computed values and the value of 6.5 meV, esti- 
mated using the two-band model together with the appropriate 
masses and screening constant, is small. 

Figure[T2]shows the exciton spectrum of InN together with 
the respective oscillator strengths for ordinary and extraor- 
dinary light polarization. Further, the contributions c^^. = 
I^CTkP °f inter-band transitions between the 3 upper- 
most valence bands v^, vg, vc and the first conduction band c 
have been analyzed. They are indicated by the relative lengths 
of the differently colored bars in the lower panel of Fig. [12] 
Obviously, as the respective contributions sum up to one, only 
the three single-particle pair states {va,c), {vb,c), and (vc,c) 
are found to contribute to the excitons in the energy interval 
shown in Fig. [12] Now, also the C exciton can be identi- 
fied. To this end, we perform a constrained calculation with a 
pair Hamiltonian including only (vc,c) transitions. The low- 
est eigenvalue of this Hamiltonian is found 2.3 meV below 
the vc-CBM gap (cf. Fig.[T2Ti. Indeed, the unconstrained pair 
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FIG. 12: (Color online) Exciton spectram of InN. The pair ener- 
gies are given with respect to the lowest gap energy. The oscillator 
strengths are given in the upper panel. Black solid bars indicate the 
average oscillator strengths in ordinary polarization direction, while 
orange bars indicate the extraordinary polarization direction. The 
cyan bars indicate eigenstates of a pair Hamiltonian, which includes 
only contributions of the and the first conduction band. The red 
dots indicate the average oscillator strength. In the lower panel the 
contributions resulting from transitions between va,vb,vc and the 
first conduction band c are analyzed (see text). 



Hamiltonian gives rise to an eigenvalue at an only marginally 
higher energy, 2.1 meV below the vc-CBM gap. The respec- 
tive exciton, however, is found to be a nearly half-half mix- 
ture of {vb,c) and (vc,c) contributions. Further, it is clear 
from Fig. [12] that the latter exciton is neither the strongest 
bound state with contributions from (vc,c) transitions, nor it 
is the first exciton visible in extra-ordinary polarization, nor 
it has the strongest oscillator strength in the latter polariza- 
tion direction. Especially the latter two points are notable, 
since they contradict the initial expectations and may affect 
the experimental assignment of excitons. The reason for the 
observed behavior is found in the k-point dependence of the 
single-particle optical oscillator strengths o= |(<Pek \Pj \ (Pvk)?'- 
While the optical transitions (v^,c) and (vb,c) at F are dipole- 
forbidden in extraordinary polarization, they are allowed and 
of comparable strength with the (vc,c) transitions at off-F k 
points. 

Finally, we address the influence of the LDAh-U approxi- 
mation on the present results. The electron and vc effective 
masses at F increase with increasing U, but underestimate 
the values calculated upon the HSE03 functional. Larger val- 
ues than U=3 eV, however, yield a vanishing or even inverted 
crystal-field splitting - contradicting the HSE03 results and 
experimental findings. Due to the accompanying crossing and 
mixing of the v^, vb, and vc bands, values higher than U=3 eV 
are not meaningful starting points. The same holds for lower 
values of U, which predict much too small gaps and electron 
masses. The effect of the underestimated electron mass can 
be estimated within the two-band model, where the binding 



energy is found to depend linearly on the reduced effective 
mass. Therefore, the computed binding energies of the A and 
B excitons may underestimate the actual values by about 30%. 

Apart from the influence of the U parameter the calcu- 
lated exciton binding energies are influenced by the screen- 
ing. Especially dynamical screening ll56il57ll may reduce the 
exciton binding. The characteristic optical phonon energies 
(HcOlo = 86 meV |54]) are much larger than R^x, so that the 
lattice polarization can fully contribute to the screening of the 
electron-hole attraction. For InN the static dielectric constant 
amounts to = 13 |75]. Using Ej for the screening instead 
of Coo would significantly reduce the exciton binding energies 
to values below 2 meV. Consequently, Wannier-Mott excitons 
cannot be observed in the currently available InN samples, 
due to their high density of free carriers which is above the 
Mott-transition energy ll6Tll62ll63ll . 



V. SUMMARY AND CONCLUSIONS 

We have introduced a new and highly efficient numerical 
approach for the solution of the homogeneous BSE for bound 
electron-hole pair excitations in non-metals. The use of an 
iterative diagonalization scheme is demonstrated to diminish 
the computational costs for the calculation of bound electron- 
hole pair states significantly, aflowing for the systematic study 
of excitonic Hamiltonians with ranks up to several hundred 
thousands. This is due to its, in comparison to the 0{N^) 
scaling of direct diagonalization methods, favorable G{N^) 
scaling. Further, a hybrid k-space sampling scheme is pre- 
sented which allows a refined sampling in selected parts of 
the BZ and simultaneously avoids an artificial localization of 
the exciton wave functions. The corresponding discretization 
of the pair Hamiltonian EVP is performed and found to yield 
a generalized EVP with a diagonal overlap matrix. 

As a first test and paradigmatic example the developed nu- 
merical approaches are applied to the Wannier-Mott exciton, 
which is solved numerically in the reciprocal space. Using 
realistic model parameters for the reduced effective mass and 
screening constant, the numerical convergence is systemati- 
cally studied with special focus on the convergence with re- 
spect to the k-point sampling. The well known analytical so- 
lution of a hydrogen-like exciton series is reproduced for the 
lowest bound pair states corresponding to the principle quan- 
tum numbers n = 1 ... 3, including the expected degeneracies 
of exciton states and the oscillator strengths. 

Further, the developed approaches are applied to the non- 
metallic materials MgO and InN, which both show Wannier- 
Mott-like excitons. The respective excitons, however, differ 
drastically in terms of their localization and consequently con- 
vergence, due to the very different electronic structures found 
for both materials. Nevertheless, the convergence trends iden- 
tified for the Wannier-Mott model are found to hold and allow 
for the extrapolation of k-point converged binding energies for 
the strongest bound excitons. If only the electronic screening 
is taken into account, the converged binding energies clearly 
exceed the experimental results. This is attributed to the influ- 
ence of the lattice polarization on the screening and therewith 
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the exciton binding. 

Further, the achieved level of convergence, together with 
the account for the complete electronic structure, allows for 
the discussion of effects, which are not addressable in k • p the- 
ory. For InN, for instance, the often neglected k-dependence 
of the single-particle oscillator strengths is shown to influence 
the polarization dependence of the lowest excitons drastically. 
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